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Abstract 

We present a systematic hierarchy of approximations for local exact-decoupHng of four-component 
quantum chemical Hamiltonians based on the Dirac equation. Our ansatz reaches beyond the trivial 
local approximation that is based on a unitary transformation of only the atomic block-diagonal part of 
the Hamiltonian. Systematically, off-diagonal Hamiltonian matrix blocks can be subjected to a unitary 
transformation to yield relativistically corrected matrix elements. The full hierarchy is investigated 
with respect to the accuracy reached for the electronic energy and for selected molecular properties 
on a balanced test molecule set that comprises molecules with heavy elements in different bonding 
situations. Our atomic (local) assembly of the unitary exact-decoupling transformation — called local 
approximation to the unitary decoupling transformation (DLU) — provides an excellent local approx- 
imation for any relativistic exact-decoupling approach. Its order- A^^ scaling can be further reduced to 
linear scaling by employing a neighboring-atomic-blocks approximation. Therefore, DLU is an efficient 
relativistic method well suited for relativistic calculations on large molecules. If a large molecule con- 
tains many light atoms (typically hydrogen atoms), the computational costs can be further reduced 
by employing a well-defined nonrelativistic approximation for these light atoms without significant loss 
of accuracy. We also demonstrate that the standard and straightforward transformation of only the 
atomic block-diagonal entries in the Hamiltonian — denoted DLH in this paper — introduces an er- 
ror that is on the order of the error of second-order Douglas-Kroll-Hess (i.e., DKII2) when compared 
with exact-decoupling results. Hence, the local DLH approximation would be pointless in an exact- 
decoupling framework, but can be efficiently employed in combination with the fast to evaluate DKH2 
Hamiltonian in order to speed up calculations for which ultimate accuracy is not the major concern. 
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I. INTRODUCTION 



The consideration of relativistic effects is essential to the proper understanding of the chem- 
istry of heavy-element containing molecules l|, |2|. The Dirac equation provides the relativistic 
quantum mechanical description of a single electron in the presence of external electromagnetic 
potentials. As the Dirac operator comprises four-dimensional operators, its eigenfunctions pos- 
sess four components. The Coulomb interaction turned out to be sufficiently accurate for the 
description of chemical phenomena and the corresponding many-electron Hamiltonian is called 
the Dirac-Coulomb Hamiltonian. 

A relativistic method based on the Dirac-Coulomb Hamiltonian is called a four-component 
method. It yields negative-energy solutions, which are pathologic and of no use in chemistry. 
Also, a large number of basis functions is required to properly describe the negative-energy 
states. These are two drawbacks of four- component methods that motivated the development 
of two-component methods, which remove the negative energy solutions and can, in principle, 
exactly reproduce the results of four-component calculations [2|. 

Several relativistic two-component methods were developed in the past decades. One of 
the widely used approaches is the second-order Douglas-KroU-Hess method (DKH2) [s], 0]. 
It employs the free-particle Foldy-Wouthuysen (FW) transformation as well as sequential 
Douglas-Kroll 6] transformations to decouple the four-component operator. Higher-order [3- 



[lol | and even arbitrary-order 
regular approximation (ZORA) 



17 



_[l6|_DKH methods have been developed. The zeroth-order 



3 i 



is another highly successful relativistic two-component 



method. Within the ZORA framework, it is particularly easy to implement the calculation of 
molecular properties; see Refs. 



20 



2l| for very recent examples and Ref. 



221 for a review. 



In recent years it was shown that, if the goal is to exactly decouple the four-component 
Hamiltonian matrix instead of the Hamiltonian operator, the formulae of exact decoupling can 
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26| method aims at exact decoupling 



be directly obtained. The Barysz-Sadlej-Snijders (BSS) 
of the free-particle FW-transformed four- component Hamiltonian matrix by a matrix operator 
of the form derived in Ref. 27|. The key matrix operator, which is used to construct the de- 



coupling transformation matrix, is obtained by solving an iterative equation. However, invoking 
the free-particle FW transformation turns out to be not necessary for the construction of the 
exact decoupling transformation in a one- step protocol. The pioneering work of this one-step 



transformation was provided by Dyall 



28l-l32| in the form of the so-called normalized elimina- 



tion of the small component (NESC) approach. Later it was generalized to the so-called exact 



two-component (X2C) method by several groups 33l-l44j. Since the DKH method also employs 
matrix techniques to evaluate the Douglas-Kroll transformation, it is also able to exactly decou- 
ple the f our- comp on ent Hamiltonian matrix. This has been shown within the arbitrary-order 

-53. 
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approach |11|-|14|. Il6l. I45l]. For reviews of these developments see Refs. 

Even within a scalar relativistic approximation, where the spin-dependent operators are ne- 
glected, the construction of the relativistic Hamiltonian matrix dominates the step of the one- 
electron integral evaluation. This is the case because the relativistic transformation is done 
within an un-contracted basis set and involves many matrix manipulations such as multiplica- 
tion and diagonalization. For a large molecule including only one or a few heavy atoms, this 
procedure would waste computational resources since the relativistic effects are highly local and 
only significant for heavy atoms. It is therefore clear that a local relativistic method is desirable. 
Even for a molecule containing only heavy atoms, such a local relativistic approach would be 
preferable as we shall see. 

The relativistic Hamiltonian operator is universal in coordinate space and a local approach 
is unfeasible at the operator level. However, since most quantum chemical calculations apply 



a linear combination of atom-centered basis functions, we may employ these basis 



unctions to 



53| . For those 



construct atomic projectors as it is done in charge and spin population analysis 
Hamiltonian matrix blocks for which relativistic corrections are important, we shall derive a 
local relativistic Hamiltonian matrix to evaluate it. The locality is then exploited in the basis- 
function space instead of coordinate space. A crucial aspect to determine will be which block 
of the Hamiltonian matrix requires a relativistic description and how to evaluate it. Undoubt- 
edly, the heavy atom diagonal blocks will require a relativistic description, while the relativistic 
description of heavy atom off-diagonal blocks depends on their contribution to physical observ- 
ables and on the feasibility of applying the relativistic description. For relativistic transformation 
techniques, such as DKH and X2C, it is impossible to apply the relativistic transformation to an 
off-diagonal block alone without information about other blocks. It requires the information of 
a full square matrix to construct the relativistic transformation. Such a difficulty disappears if 
one employs an operator-based relativistic method such as ZORA. The off-diagonal block of the 
Hamiltonian matrix can be directly computed from relativistic integrals. However, the ZORA 
method provides a poor description of atomic inner-shell orbitals, which carry the largest part 
of the relativistic effect. In this sense, the exact decoupling methods are the best candidate for 
applying a local relativistic scheme since they do not neglect important relativistic effects. 
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All current exact decoupling methods employ the transformation technique. It is therefore 
not trivial to apply the relativistic description to heavy-atom off-diagonal Hamiltonian matrix 
elements. However, in this work we propose that if the local approximation is applied to the 
decoupling transformation instead of the Hamiltonian matrix, one could obtain results which 
include the relativistic description to heavy atom off-diagonal blocks. We also note that the 
main effort for the construction of the relativistic Hamiltonian matrix with exact-decoupling 
methods is due to the evaluation of the decoupling transformation matrices. Therefore, a local 
approximation to the decoupling transformation would also lead to a significant reduction of the 
computational cost. 

II. RELATIVISTICALLY LOCAL STRUCTURE OF THE HAMILTONIAN MATRIX 

The following notation will be used throughout this paper. Upper case labels A, B, C denote 
heavy atoms which require a relativistic description, while lower case labels a, 6, c denote light 
atoms for which a (standard) nonrelativistic description can be assumed to be accurate enough 
based on numerical evidence compiled within computational chemistry in the past decades. To 
be more precise, the former labels refer to nuclei with high nuclear charge numbers, while the 
latter refer to those with low nuclear charge numbers (typically with nuclear charges of less than 
two dozen protons). Of course, the distinction whether a nucleus is considered heavy or light 
bears some sort of arbitrariness and it will be made on the basis of the accuracy required in 
a quantum chemical calculation. In general, there are then five different types of Hamiltonian 
matrix blocks Haa, Hab, Hau, Haa, and Hab that need to be considered. 

The evaluation of the relativistic Hamiltonian matrix with exact-decoupling approaches re- 
quires a set of square matrices. They are the nonrelativistic kinetic energy matrix T, 




^(A,|p2|Aj), 




the external potential matrix V, 



V.J = (A.|V|A,), 



(2) 



and the overlap matrix 5", 




(3) 



as well as an additional relativistic matrix W 




{\i\(T -pVa -plXj), 



(4) 



5 



(note that all expressions are given in Hartree atomic units, in which the rest mass of the 
electron takes a numerical value of one). In the above equations, Aj denotes the i-th 2-spinor 
atom-centered basis function, V the external potential, c the speed of light, cr the vector of 
Pauli spin matrices. The basis functions Aj may be grouped according to the atomic nucleus 
to which they are assigned, which is key to the local approaches discussed in the following. 
The matrix representation of the relativistic Hamiltonian is then evaluated as a function of the 
above-mentioned matrices 

H = H{T,V,W,S). (5) 

A straightforward idea for a trivial local approximation is to apply the evaluation of the 
relativistic Hamiltonian matrix only to atomic (i.e., diagonal) blocks 

Haa = Haa^Taa, Vaa, Waa, Saa), (6) 

(instead of to the full matrix, but where A may also represent a group of atoms) and to ignore 
the relativistic correction to all off-diagonal blocks 

Hab = Tab + Vab- (7) 
Adding more atoms to a group A will improve the accuracy, but the computational advances 



will be lost when a group A becomes very large. This approach was employed in Ref. 
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55| 



and explored in detail in Ref. 56| for the low-order DKH method. We refer to it as the diagonal 
local approximation to the Hamiltonian (DLH) matrix. It is clear that the DLH approximation 
can also be applied to relativistic exact-decoupling approaches. The DLH approximation works 
well at large interatomic distances, but the difference to a reference energy increases with shorter 
distance. 

However, relativistic corrections to off-diagonal blocks are also important especially when 
inter-atomic distances are short. The necessity of applying the relativistic description to off- 



diagonal 
In Ref. 



ocks has already been discussed in Ref. 



571 ] with a so-called two-center approximation. 



57|, the DKH transformation was applied to pairs of atoms {A © B) 

Ha®B = H{TA®B-,yA®B-,WA®B-,SA®B)-, 

with 

I, (9) 

\ Hba Hbb 
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to obtain a relativistically corrected off-diagonal Hamiltonian matrix Hab- Unfortunately, such a 
two-center approximation introduces an inconsistent treatment to diagonal blocks. For instance, 
Haa can be obtained from either the pair {A © B) or from another pair {A(B C). 

So far, we have only discussed the AA and AB blocks which, by definition, require a relativistic 
description for both atoms. Since relativistic corrections for light atoms are very small and may 
be neglected, Haa and Hab can be approximated by the nonrelativistic Hamiltonian matrices 

Haa = Taa + Vaa, (10) 
Hab = Tab + Vab. (11) 

A relativistic scheme for the heavy-light hybrid off-diagonal blocks HAa should, however, be 
considered explicitly. 



III. LOCAL DECOMPOSITION OF THE X-OPERATOR 



The general expression of exact decoupling transformations can be written as 



U 



( 



\ 



-X z- 







ULS\ 






jjSS j 



(12) 



Z~^ and Z~ are two-component unitary operators, and only Z~^ is required for the evaluation 
of the electrons-only Hamiltonian. The X-operator generates the electronic small- component 
functions ip'^ from the large-component functions ip~l 



(13) 



Two-dimensional Hamiltonians are then obtained by applying the unitary transformation of Eq. 
(fT2|) to blockdiagonalize ('6(i') the four-dimensional Dirac-based Hamiltonian D 



H'"^ = UDUl 



(14) 



It yields both the electrons-only Hamiltonian H (as the upper left block in H'"^) and the one for 
negative-energy solutions, although the latter will be discarded. 



As discussed by Dyall in Ref. 



29| for the NESC method, the Breit-Pauli approximation to 



the X matrix, which is used to construct the relativistic transformation matrix, is 

X = L 



(15) 
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where / denotes the identity matrix. Such an approximation is not variationally stable and thus 
cannot be used for variational calculations. However, one could employ this approximation for 
light atoms only 



(16) 



since the corrections to light atoms are small and may not affect variational calculations. This 
idea was suggested in Ref. 42] for the X2C method but no results were presented. 

Although the lowest level approximation to X Eq. (|T5|) did not provide useful results, the 
atomic approximation to the X matrix gave results with very small errors. As discussed by 



Dyall in Ref. 



30 



3l|, the X matrix can be approximated as the direct sum of atomic blocks 

X = Xaa®Xbb®--- (17) 



This approximation resembles the linear combination of atomic four-spinors ansatz in four 
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in the sense of an 



component calculations. The same approach was also employed in Ref. 
atomic approximation to the projection on electronic states. If the electronic states are further 
transformed to a two-component picture by a renormalization matrix, this gives rise to a local 



X2C method. As discussed in Refs. 



38 



42| . the local approximation to the X matrix Eq. ([T 



works well for spectroscopic constants of diatomic molecules. 

However, there are some drawbacks of the local-X matrix approximation to the X2C method. 
Firstly, it only reduces the cost for the evaluation of X matrix, while the exact- decoupling 
transformation matrices are still of molecular dimension as the renormalization matrix is not 
atomic block diagonal even if the X matrix is. The computational demands for the evaluation of 
the relativistic approximation are thus still tremendous for large molecules. Secondly, the X = I 
approximation did not provide a satisfactory treatment of the Aa blocks, since it may suffer from 
variational collapse and Eqs. fIlOp and fill I) cannot be fulfilled within this approximation. 



IV. LOCAL APPROXIMATIONS TO THE EXACT-DECOUPLING TRANSFORMA- 
TION 

In this article, we suggest an local approximation to the exact-decoupling transformation. We 
may approximate the unitary transformation by taking only the 'atomic' diagonal blocks (all 
off-diagonal blocks are then set to zero) 

U = Uaa®Ubb®--- , (18) 
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where the atomic unitary transformations Uaa are obtained from the diagonahzation of the 
corresponding Daa blocks of matrix operator D. We denote this approximation to U as the 
diagonal local approximation to the unitary decoupling transformation (DLU). 

If we now substitute the approximate unitary transformation of Eq. f ll8p into Eq. f ll4p . the 
diagonal blocks, 

H'^\ = UaaDaaU\a^ (19) 
turn out to be the same as in the DLH approach, while the off-diagonal blocks then read 

i/^lj = UaaDabU^bb^ (20) 

which is to be compared with the expression 

H% = J2^AiDjjUl^, (21) 
i,j 

where U has not been approximated (here, / and J run over all atomic blocks). Hence, the 
DLU approach also introduces a relativistic description to off-diagonal 'interatomic' blocks when 
compared with the DLH approximation. 

The cost for the assembly of the unitary transformation U within the DLU approach is then of 
order A^, where measures the system size, namely the number of atoms. It is linear scaling since 
the atomic approximation is directly applied to the unitary transformation. However, the next 
step of applying the unitary decoupling transformation to obtain the relativistic Hamiltonian 
matrices is no longer linear scaling. If no local approximation was applied for the unitary 
decoupling transformation, according to Eq. fl2T]) where the summation indices /, J run over all 
atomic blocks, the calculation of Hab matrix will require 2N'^ matrix multiplications. Since the 
number of matrices to be calculated is of order A^^, the total cost of the relativistic transformation 
without local approximation is therefore of order A^^. If the DLU approximation is applied, no 
summation will be needed and only two matrix multiplications will be required for each heavy- 
heavy block Hab- The total cost is then of order A^^. 

If the distance of two atoms A and B is sufficiently large, the relativistic description of Hab 
can be neglected. Thus, we may define neighboring atomic pairs according to their distances. 
Then, only the Hamiltonian matrix H'^^ of neighboring pairs requires the transformation 

H% = UaaDabUIj,, VA, 5 being neighbors, (22) 
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whereas all other pairs are simply taken in their nonrelativistic form. Since the number of 
neighboring pairs is usually a linear function of system size. The total cost is then reduced to 
order A^. 

To investigate the relativistic correction to hybrid H^a blocks of the electrons-only Hamilto- 
nian H, we need to consider the explicit structure of the unitary decoupling transformation, 



T (W-T) j \ t/^^'t 



(23) 



where only the electronic part of the exact-decoupling transformation represented by U^^ and 
U^^ [which are the upper part of Eq. (fT2!) ] is employed. From Eq. (123!) it is easy to see that the 
electrons-only Hamiltonian reads 

— JjLSrpjjLL,\ _|_ jjLLrpjjLS,\ _ jjLSrpjjLS,\ 

+U^^VU^^'^ + U^^WU^'^'l (24) 
In the nonrelativistic (NR) limit, where the speed of light approaches infinity, we have 



If we insert them into the relativistic electrons-only Hamiltonian fl24j) . we arrive at the nonrela- 
tivistic Hamiltonian matrix, 

^NR = 2- + y_ (26) 

This suggests a solution for the relativistic treatment of the hybrid heavy-light blocks Hau, 
because we may set all light-atom diagonal blocks in the unitary transformation to identity 
matrices 

[/i;^ = / and t/f/ = /. (27) 

To reproduce the nonrelativistic Hamiltonian matrices of light-atom-only blocks, we must set 
the W matrix of light atoms to zero, 

Waa = 0. (28) 

The WAa matrices are also set to zero 



WAa = 0. 
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(29) 



We denote this approximation for the treatment of hybrid blocks as the DLU(NR) approach. 
If we do not introduce any approximations to the W matrices, the transformation will yield a 
Hamiltonian that differs from the nonrelativistic Hamiltonian matrix on the order of which 
corresponds to the order achieved in the Breit-Pauli approximation. We therefore denote it as 
the DLU(BP) approach. 

The final local approximation to the unitary matrices is then 

= f/ii©f/|£©---f/i^©f/^,^©--- , (30) 

U^s ^ ^Ls^ © f/|| © . . .f/f/ © Ut,' © ■ ■ ■ , (31) 

Where the atomic unitary transformation matrices of heavy atoms are evaluated by the rela- 
tivistic exact decoupling approaches. The explicit expressions of Hamiltonian matrix blocks H^a 
and Hab are 

Haa = UaaTaaUaa'^ + UaaTaaUaa^ 

JjLSrp jtLS,] . ttLLj/ ttLL,-\ 

-^AA-I-AA^AA + '-'aA^AAUaA 

+f/iiW^^At/^f , (32) 

Hab = UaaTabUI^b'' + ^aaTabUbb'' 
-UaaTabU^b^ + UaaVabU^b^ 

+U^'aWabU^'b\ (33) 

If the DLU(NR) approach is employed, the unitary transformation matrices of light atoms are 
simply set to identity matrices and the explicit form of Hau, Haa, and Hab is 

HAa = UAA^Aa + UAA^Aa, (34) 
Haa = Taa + Ka, (35) 
Hab = Tab + Vab. (36) 

In the DLU(BP) approach, however, they read 

HAa = U^ATAa + U^AVAa + U^^WAa, (37) 
Haa = Taa + Vaa + W^aa, (38) 
Hab = Tab + Vab + Wab- (39) 
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The relativistic transformation of a property operator P 13|, ll6|, |45|, |59|, |60| requires an 
additional matrix Q and its matrix elements are given by 



The relativistic corrected property matrix then reads 

P = IjLLpjjLL,] ^ jjLSQjjLS,t_ 



(40) 



(41) 



V. DLU EVALUATION WITHIN EXACT-DECOUPLING APPROACHES 



A. The X2C approach 



The X2C method employs 



U- 



X2C 



-X 



(42) 



as the exact-decoupling unitary transformation. The unitary operators for the electronic states 
are then 

1 



TjLL 



Vi + xtx' 

X2C ~ '-^X2C" 



--X\ 



(43) 
(44) 



For an actual implementation, we need their matrix representations in a finite non-orthogonal 
space of the atom-centered basis functions. The X-matrix is evaluated by diagonalizing the 
following generalized eigenvalue equation 




S 



2^^ 




(45) 



V T 

^ T {W-T) , y^s I \ " 2^W \-5 

where C\ and denote the large- and small-component molecular-orbital coefficients, respec- 
tively. We obtain the X matrix by the following equation 

X = Ct{Ciy\ (46) 

The matrix representation of and U^^ is then given by 

^1/2(^-1/2^^-1/2^-1/2^-1/2^ (47) 

(48) 



LL 
X2C 

LS 



'^X2C ~ '^X2C^ ; 
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with S defined as 

S = S + ^X^TX. (49) 

In our atomic approximation to the unitary transformation U, we need to evaluate the above 
equation for each 'atomic' block of the original Hamiltonian in order to obtain the diagonal 
blocks f/_4^ and U^j^. The decoupling transformations are then assembled according to Eqs. (130!) 
and flHTl) and the relativistic Hamiltonian matrix follows from Eq. f l2^ . 



B. The BSS approach 

In the BSS approach, all matrices are transformed to an orthonormal basis-function space first. 
It is convenient to calculate the transformation matrix K by diagonalizing the nonrelativistic 
kinetic energy matrix 

TK = SKt, (50) 

since the eigenvalues t can be used for later evaluation of the free-particle FW transformation. 
The eigenvector matrix K has the following properties: 

K^SK = I and K^TK = t. (51) 

The free-particle FW (fpFW) transformation features four diagonal block matrices 

7-rLL,fpFW TTLS,{pFW \ I , I Ep+c^ I Eq-c^ 

t/o= I = 1 (52) 




with Eq = a/ 2tc^ + and / = jt. It is then applied to yield a transformed four-dimensional 
Hamiltonian matrix Dq 

[ K oV [ V T \ ( K , ^ 

Do = Uo\ \Ul (53) 

\o K J y T (w-T) J yo K 

Next, the exact-decoupling BSS transformation. 



-R 



(54) 
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which has the same structure as the exact-decouphng transformation f/x2C; is apphed. The R 
matrix is obtained by diagonahzing the free-particle FW transformed four-component Hamilto- 
nian matrix Dq and employing a similar relation to the one in Eq. (H6|) . 

After the exact-decoupling BSS transformation has been carried out, the Hamiltonian matrix 
is back-transformed to the original (non-orthogonal) basis 

H=\ U.DoUl . (55) 

Consequently, we obtain the matrix form of U^^s ^bss BSS approach as 

[/llg = (J^-l)t(l + i?i?t)-|(f/WpFW^^tf;SL,fpFW^^t^ (56) 
f/||g = (K-l)t(l + i?i?t)-|(f/WpFW^^t^55,fpFW)^t. (57) 

The expressions for the relativistic Hamiltonian matrix (and property matrices) are then as given 
above. 



C. The DKH approach 

The DKH approach features the same initial transformation as the BSS approach to obtain 
the transformed Hamiltonian matrix Dq. Subsequent decoupling transformations are expressed 
as 

1 

t/M =11 Uk (58) 



with the generalized parametrization of the Uk [9| 

[m/k] 

Uk=Yl (59) 

i=0 

where are anti-hermitian matrix operators and m is related to the order of the DKH expan- 
sion. Since the expression of is too lengthy to be presented here, we refer the reader to Ref. 



16| for details. If m is large (strictly, if it approaches infinity), exact decoupling is achieved. 
Usually, a low value for m is often sufficient for calculations of relative energies and valence-shell 
properties (i.e., m=l for the original DKH2 approach). 

With the complete decoupling DKH transformation written as 

= I rjSL rjSS ' " ^^^^ 
'^DKH "^DKH 
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the matrix forms of U^^^n f/oxH 

U^^^ = (i^-l)t(f/LL,M^LL,fpFW^^LS,M^5L,fpFW)^t^ (gl) 
t^DKH = (A-l)t(f/^^'Mf/^^'fpFW^f;L5,Mf;5S,fpFW)^t_ ^g2) 

We should note that the finite-order DKH Hamiltonian matrices are not the resuh of directly 
applying the decoupling transformation. To give consistent results, the transformation in Eq. 
fl23l) is only used to obtain the off-diagonal blocks while the diagonal blocks are replaced by the 
traditional finite-order DKH Hamiltonian. However, for high-order calculations, the differences 
of with/without replacing the diagonal blocks are very small. 



VI. RESULTS AND DISCUSSION 



As we have shown, there exists a systematic hierarchy of several levels for the local relativistic 
approximation to Hamiltonian matrix elements. In this section, we study the accuracy of the 
different levels on a test molecule set and introduce the following notation to distinguish the 
results: 

1. FULL: Full molecular relativistic transformation, no approximation. 

2. DLU: Diagonal local approximation to decoupling transformations, U^^ and U^^ are 
evaluated for each atomic (diagonal) block. However, three variants are possible: 

i) DLU(ALL): Calculate relativistic transformation for all atoms (all atoms are treated 
as 'heavy' elements). I.e., apply the 'atomic' unitary transformation to the diagonal 
and to all off-diagonal blocks of the Hamiltonian 

ii) DLU(NB): Calculate relativistic transformation for all atoms (all atoms are treated 
as 'heavy' elements), but apply the relativistic transformation only to blocks of the 
Hamiltonian assigned to neighboring atoms (and to the diagonal blocks, of course) to 
achieve linear scaling. 

iii) DLU {ABC, BP): Distinguish heavy and light atoms, where ABC denote the heavy 
atoms and apply the BP approximation to the light atoms. 

iv) DLU(yli?C,NR): Distinguish heavy and light atoms, where ABC denote the heavy 
atoms and apply the NR approximation to the light atoms. 
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3. DLH: Diagonal local approximation to Hamiltonian matrix blocks. The relativistic cor- 
rection for off-diagonal 'interatomic' blocks are neglected. 

i) DLH(ALL) : Calculate a relativistic description for all diagonal 'atomic' blocks. 

ii) DLH(Ai?C,NR/BP) : Only heavy-atom blocks are considered for the relativistic 
description. 

Note that this classification system does not account for all possible local relativistic approxima- 
tions: i) A could also represent a group of atoms and ii) additional off-diagonal approximations 
may also be included. 

We have implemented the local relativistic schemes discussed so far for the X2C, BSS, and 
DKH exact- decoupling approaches into the MoLCAS programme package 6l|, into which we 
have recently implemented X2C, BSS, and the polynomial-cost DKH schemes 52|. In our 
scalar-relativistic calculations the definition of the W matrix was 

W^j = ^{\^\p■Vp\\J), (63) 

and Aj now refers to scalar basis functions instead of 2-spinor basis functions. All other expres- 
sions derived in this work are the same for the scalar-relativistic X2C, BSS, and DKH variants. 
We should mention that we have not considered the transformation of the two-electron integrals 
as it is common in most exact-decoupling calculations. However, the local methods discussed 
here could be of value also for one of the approximate methods available to transform the two- 
electron integrals (note that a full-fledged transformation of the two-electron integrals would 
render the application of exact-decoupling methods inefficient and a four-component approach 
should be employed instead). 

For our test molecule set, we should rely on closed-shell molecules with heavy elements in 
the vicinity of other heavy elements and of light elements. Four molecules were selected to test 
the validity of the different local schemes: I^, WH6(PH3)3, W(CH3)6 and Pbg"*". Moreover, two 
reactions 

1+ ^ 1+ + h (64) 

WH6(PH3)3 ^ WH6(PH3)2 + PH3 (65) 

were chosen to study reaction energies. Their structures (see Fig. [1]) have been optimized with 



the TURBOMOLE program package 



62l | employing the BP86 density functional 63|, |6J] with a 
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valence triple-zeta basis set with polarization functions on all atoms 65|] in combination with 
Stuttgart effective core potentials for W, Pb, and I as implemented in Turbomole. 






3.00 



I 

12.70 
'I 




;.70 
'I 



3.00 




FIG. 1: Structures of the test molecule set (white spheres are hydrogen atoms; selected bond lengths 
are given in A). 

All molecules are closed-shell and for our analysis all- electron Hartree-Fock calculations were 
then performed with the MoLCAS programme package. Since our aim is to test the validity of 
local relativistic approximations instead of obtaining accurate results which can be compared 
to experiments, we do not need to consider electron-correlation effects in our calculations. The 



66 



-68 



were employed at a double-zeta 



all-electron atomic natural orbital (ANO) basis sets 
level for all atoms. The exact-decoupling DKH calculations were performed at 35-th order, since 
higher order results did not improve the accuracy for our purpose so that this high DKH order 
can be considered exact. Results of finite-order (standard) DKH2 calculations are also given for 
the sake of comparison. 

The computation times on the Intel 17-870 CPU for the different local schemes are compared 
at the example of Pbg"*". Table [T] presents the timings for the evaluation of the relativistic one- 
electron Hamiltonian matrix operator, (the picture change transformation of property integrals 
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TABLE I: Computation times (in seconds) for the relativistic one-electron Hamiltonian in various local 
schemes (for Pbg"*"). 

DKH35 BSS X2C DKH2 

DLH(ALL) 1216 6 6 1 

DLU(ALL) 1225 15 15 11 

FULL 124900 564 520 92 



are not considered in the measurement of the CPU time). For comparison, the computation times 
for the one-electron and for the electron repulsion integrals are 8s and 29744s, respectively. They 
are, of course, the same for the different relativistic approaches. Since the relativistic transforma- 
tions are performed on each atomic block in the DLH scheme, its computational costs are much 
smaller than the FULL scheme. The ratio of DLH to FULL is almost a constant (around 1% in 
this case) for all relativistic approaches. The DLU scheme has an extra transformation step for 
off-diagonal blocks in addition to the DLH scheme, it therefore requires a constant time (about 
10 second in this case) more than the DLH scheme, which is still negligible when compared to 
the FULL scheme. The computation times for the evaluation of both one-electron integrals and 
relativistic transformations (except for DKH35) are smaller than that for two-electron integrals. 
The relativistic transformation is usually not the bottleneck step. However, if density functional 
theory calculations — with a pure generalized-gradient-approximation functional in combination 
with some acceleration scheme like density fitting or fast multipole approximations — instead of 
Hartree-Fock calculations were performed, the computation time of the relativistic transforma- 
tion would be significant and the local approximation would then be required. Of course, this 
also holds especially for much larger molecules with light atoms (see also discussion below). For 
the high-order DKH method, like DKH35 employed here, the computational cost of the FULL 
scheme even exceed that of the two-electron integrals and therefore a local approximation is 
mandatory. 

We performed DLU(ALL) and DLH(ALL) calculations for all four molecules in order to 
analyze the general error that is introduced by these local approximations. The differences 
of total electronic energies with respect to the FULL reference calculation without any local 
approximation are presented in Table [Tll Since the errors of total energies are rather small, 
all values are given in milli-Hartree (niH) atomic units. The errors in total energies can be 
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understood as errors on atomization energies that one would observe when atomization energies 
are to be calculated with one of the local approaches compared to the FULL one. 



TABLE II: Total energy differences introduced by local approximations DLU(ALL) and DLH(ALL) 
with respect to the FULL reference energy. All values are in milli-Hartree atomic units. 



Method Approximation 



I 



WH6(PH3)3 



W(CH 



3)6 



DKH35 FULL 

DLU(ALL)-FULL 
DLH(ALL)-FULL 



-35562819.0554 
-0.0119 
0.2471 



-17170991.3602 
-0.0126 
-30.8580 



-16375243.0202 
-0.0204 
-6.3823 



-187905906.7440 
0.0376 
-2.2619 



BSS FULL 

DLU(ALL)-FULL 
DLH(ALL)-FULL 



-35562831.5438 
-0.0119 
0.2461 



-17171044.3875 
-0.0125 
-30.8303 



-16375296.0488 
-0.0204 
-6.3749 



-187907231.0132 
0.0387 
-2.2477 



X2C FULL 

DLU(ALL)-FULL 
DLH(ALL)-FULL 



-35561841.0814 
-0.0133 
0.2629 



-17171053.6272 
-0.0306 
-33.3336 



-16375311.3620 
-0.0367 
-7.0463 



-187911251.8227 
0.0523 
-2.8362 



DKH2 FULL 

DLU(ALL)-FULL 
DLH(ALL)-FULL 



-35553068.6409 
-0.0083 
0.2897 



-17158776.3373 
0.0011 
-33.4347 



-16363031.9017 
0.0143 
-6.9555 



-187712096.3589 
0.0122 
-3.3895 



For the I5 molecule, which possesses a pseudo-one-dimensional structure with iodine atoms 
connected to at most two other iodine atoms, the DLU(ALL) scheme produces an error of about 
0.012 mH for both DKH and BSS, while X2C has a sightly larger error of 0.013 niH. Interestingly, 
it is smallest for DKH2, i.e., only 0.008 mH. The DLH(ALL) local scheme produces relatively 
large errors on the total electronic energies, they are 0.24 mH for DKH and BSS, 0.26 mH 
for X2C, and 0.29 mH for DKH2. Clearly, DLU(ALL) provides more accurate total energies 
than DLH(ALL), as it includes the relativistic form of the atom-other-atom off-diagonal blocks, 
while DLH(ALL) only takes into account the atom-same-atom diagonal blocks. The errors 
of DLH(ALL) are about 20 times larger than those of DLU(ALL) for the I^ molecule. This 
indicate that the relativistic corrections to off-diagonal blocks are quite important. Apparently, 
DLU(ALL) provides an accurate scheme for taking the off-diagonal relativistic description into 
account. 
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The Pbg"*" molecule features a more compact structure where more than two heavy atoms 
bind to a central lead atom. Moreover, lead atoms are heavier than iodine atoms and relativistic 
effects are more pronounced. As we can see from Table [TTl errors of Pbg"*" are all larger than 
those observed for the molecule. The deviations of DLU(ALL) electronic energies from the 
FULL reference energy are 0.038 mH for BSS and DKH, while X2C is still slightly larger with 
a value of 0.052 niH (again DKH2 features the smallest deviation of only 0.012 mH). However, 
deviations of DLH(ALL) results from the reference have become significantly larger, they are 
2.27, 2.25, and 2.84 mH for DKH, BSS, and X2C, respectively. Hence, the errors of DLU(ALL) 
are roughly 60 times smaller than those observed for DLH(ALL). 

Considering now molecules that also contain light atoms, the relative accuracy of DLU(ALL) 
vs. DLH(ALL) changes dramatically. First of all, we note that our DLU(ALL) scheme preserves 
its accuracy and the deviations from the FULL reference electronic energy are in between those 
obtained for If and Pbg"*" with the same trends as already discussed. The errors of DLH(ALL) 
energies obtained for the WHe(PH3)3 molecule are, however, larger than 30 mH. I.e., they are 
significantly larger than those observed for Pbg"*". For W(CH3)6, the deviations from the reference 
are around 6 mH, which are also still larger than for Pbg"^ and certainly nonnegligible. Note that 
the DLH(ALL) errors obtained for WH6(PH3)3 and W(CH3)6 are not all below those of the Pbg"*" 
molecule although the nuclear charge number of tungsten is smaller than that of lead. This is in 
contrast to the results obtained for our DLU(ALL) scheme. Since DLH neglects any relativistic 
correction to off-diagonal 'interatomic' blocks, the large errors must be due to the strong coupling 
between tungsten and its surrounding lighter atoms (phosphorus and carbon). Apparently, 
if a heavy-element-containing molecule has several strong bonds to this heavy element, the 
off-diagonal blocks of neighboring atom pairs can have significant contributions to the total 
electronic energy. Therefore, the DLH scheme may then introduce non-negligible errors. 

Since DLU(ALL) gives a balanced description for all blocks, it will still work in those cases 
where DLH(ALL) fails. DLU can therefore even be recommended for calculations of molecular 
systems in different geometries, as for example required for the calculation of potential energy 
surfaces. DLH would be unfit for this case as it will introduce large errors when the inter-atomic 
distances become short. 

DLU (ALL) also works well for the finite-order (standard) DKH2 method. The errors are even 
smaller compared to exact-decoupling approaches. However, note that the reference energies of 
the different relativistic approaches are not the same. Since the DLU (ALL) scheme only discards 
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the off-diagonal bloclcs of the decouphng transformation compare to the FULL reference scheme, 
the small errors in the total electronic energies indicate that neglecting off-diagonal terms in 
the decoupling transformation produces negligible errors for the total energies. This conclusion 
holds for both exact-decoupling and finite-order DKH approaches. With the DLH(ALL) scheme, 
however, the errors in the DKH2 calculation that are due to the DLH approximation are of the 
same order as those introduced by the finite, i.e., second order when compared to the exact- 
decoupling approaches. This observation leads to the very important conclusion that none of 
the exact-decoupling approaches is superior to a finite-order DKH approach if one employs a 
DLH scheme as this would render the accuracy gained by exact decoupling meaningless. In 
fact, if DLH is for some reason the local approximation of choice, a standard DKH2 approach 
would yield a more efficient relativistic scheme and nothing would be gained by achieving exact 
decoupling. 

Besides total energies, we also study the effects of local approximations on molecular proper- 
ties. We calculated the picture change corrected quadrupole moments. The property matrices 
are transformed according to Eq. fHT]) . The values for the components of the quadrupole mo- 
ment depend on the orientation of the molecule. Therefore, we present the spherically averaged 
quadrupole moments defined as ^J (Qxx + Qyy + Qzz)/'^^ Table IIIII Where Qxx denotes the 
XX (diagonal) component of the quadrupole moment. The spherically averaged quadrupole 
moments are isotropic and thus independent of molecular orientation. 
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TABLE III: Spherically averaged quadrupole moments for our test set of molecules. All values are given 
in Hartree atomic units. 



Method 


Approximation 




WH6(PH3)3 


W(CH3)6 




DKH35 


FULL 


83.84331 


57.61285 


57.29508 


173.68894 




DLU(ALL)-FULL 


-0.00009 


0.00018 


0.00025 


0.00416 




DLH(ALL)~FULL 


0.02176 


-0.02448 


-0.00218 


-0.20252 


BSS 


FULL 


83.84329 


57.61285 


57.29507 


173.68801 




DLU(ALL)-FULL 


-0.00009 


0.00018 


0.00025 


0.00414 




DLH(ALL)-FULL 


0.02177 


-0.02449 


-0.00218 


-0.20248 


X2C 


FULL 


83.84406 


57.61240 


57.29499 


173.68292 




DLU(ALL)~FULL 


-0.00017 


0.00029 


0.00027 


0.00590 




DLH(ALL)-FULL 


0.02199 


-0.02483 


-0.00197 


-0.21653 


DKH2 


FULL 


83.85310 


57.61393 


57.29674 


173.80677 




DLU(ALL)-FULL 


0.00009 


-0.00012 


0.00018 


-0.00003 




DLH(ALL)-FULL 


0.02196 


-0.02458 


-0.00186 


-0.20760 



The relative errors of DLU(ALL) are all below 0.001%, this upper limit can be further reduced 
to 0.0001% if we do not include Pb9+. Hence, DLU(ALL) is not only a good approximation 
for total electronic energies but also for molecular properties like quadrupole moments. The 
errors introduced by DLU(ALL) are all negligible, whereas DLH(ALL) shows large errors for 
quadrupole moments. While for total electronic energies, the largest error of the DLH(ALL) 
scheme shows up in WHe(PH3)3 and the smallest in I^, it is, for the quadrupole moments, largest 
for Pbg"*" and smallest for W(CH3)6. This different behavior of DLH(ALL) indicates again the 
important contribution of off-diagonal 'interatomic' blocks, whose importance must be different 
for electronic energies and for properties. In our DLU(ALL) scheme, I^ has the smallest error and 
Pbg"*" the largest one for both, energies and properties. Therefore, DLU(ALL) turns out to give 
a consistent description of relativistic effects, while DLH(ALL) yields an unbalanced relativistic 
treatment. We should stress that we have also investigated other properties (radial moments, 
electric field gradients, contact densities) drawing basically the same conclusions (which is the 
reason, why we do not report those data here). 
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For both, energies and quadrupole moments, in the DLU(ALL) scheme, DKH and BSS yields 
almost the same deviations from the reference and X2C always has larger ones than DKH and 
BSS. This similarity of DKH and BSS can be well understood since both of them employ the 
free-particle FW transformation as well as the subsequent exact-decoupling transformation. The 



close results of DKH and BSS were also observed in Ref. [52|. The slightly larger errors of X2C 
indicate that the off-diagonal blocks of the X2C exact-decoupling transformation have slightly 
larger contributions compared to DKH and BSS. However, such differences of errors are too small 
to be considered a serious drawback for the X2C approach; they are smaller than the differences 
of the total values. We may use WH6(PH3)3 as an example. The difference of total isotropic 
quadrupole moments (57.61285—57.61240) is 0.45 milli atomic units, while the difference of 
errors by DLU(ALL) (0.00029—0.00018) is 0.11 milli atomic units. The errors introduced by the 
DLU(ALL) approximation are smaller than the differences between exact-decoupling approaches. 
This is also the case for the differences between DKH2 and exact-decoupling approaches. As we 
can see, the difference between DKH2 and BSS of total isotropic quadrupole moments (57.61393- 
57.61285) is 1.08 milli atomic units, which are much larger than the DLU(ALL) errors. 

For chemical purposes, comparisons of total electronic energies are meaningless, because chem- 
ical reaction kinetics and thermodynamics are governed by energy differences. We therefore 
calculated the electronic reaction energies of the two reactions mentioned above. The results are 
given in Table IIVI and Table |Vl The errors of relative energies are less than those of absolute 
total energies. For example with the DKH approach, the error of is reduced from 0.012 mH 
to 0.002 mH for DLU(ALL), 0.24 mH to 0.12 mH for DLH(ALL). For the WH6(PH3)3 molecule, 
the error is decreased from 0.0126 mH to 0.0006 mH for DLU(ALL), 30.9 mH to 19.4 mH for 
DLH(ALL). This also holds for other relativistic approaches. 

TABLE IV: Electronic Hartree-Fock reaction energies for the reaction I^" — > + I2. All values are in 
milli-Hartree (note that 1 mH is equivalent to about 2.6 kJ/mol). 

DKH35 BSS X2C DKH2 

FULL 13.9797 13.9798 13.9739 13.9550 

DLU(ALL)-FULL 0.0022 0.0022 0.0026 0.0011 

DLU(NB)-FULL 0.0040 0.0040 0.0044 0.0027 

DLH(ALL)-FULL -0.1190 -0.1190 -0.1214 -0.1231 
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Next, we study the neighboring-atomic-block (NB) approximation at the example of the 
molecule. As shown in Figure [H it is composed of a chain of iodine atoms and thus suitable 
for applying the neighboring approximation. Only the pairs of atoms which are connected with 
a bond displayed in Figure [T] are considered as neighbors. The neighboring groups could also 
be determined by introducing a cut-off distance parameter. Then, the pair of atoms which has 
shorter distance than the given parameter is counted as a neighboring pair. The relativistic 
transformations are now only applied to such neighboring pairs according to Eq. (122|) within 
the DLU scheme. As we can see from Table HVl the neighboring approximation DLU(NB) gives 
very small errors for relative energies. Therefore it is a very good (additional) approximation 
within our DLU scheme. Hence, with DLU(NB) the computational cost will become linear 
scaling, while the DLU(ALL) scheme has an order-A^^ scaling, which may be a bottleneck for 
calculations on very large molecules. But be aware that the success of the NB approximation 
depends on how one defines the neighboring atoms. If the cut-off distance is too small, the DLU 
scheme will be reduced to the DLH scheme since no neighbors would be found. 

Reduction in terms of computational costs can also be achieved by employing the BP or 
NR approximation to the light atoms. It should be a good approximation if a molecule has 
many light atoms such as hydrogen. We study these approximations for the phosphine ligand 
dissociation energy of WH6(PH3)3. The results are presented in Table IVl where DLU(WP,NR) 
denotes the NR approximation applied to the H atoms while W and P atoms are relativistically 
described. By contrast, DLU(W,NR) denotes application of the relativistic transformation only 
to the W atom, while for P and H the NR approximation applies. From Table |V] we can see 
that the results of the BP approximation are very bad, especially if the BP approximation was 
employed for the phosphorus atoms. The BP approximation to P and H atoms gives errors 
which are roughly seven times the value for the reactions energy. Hence, this is approximation 
turns out to be completely useless. Such huge errors stem from the variational collapse of the 
BP approximation. 
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TABLE V: Relative energy difference of the reaction WH6(PH3)3 — > WH6(PH3)2 + PH3. All values 



are in milli-Hartree. 





DKH35 


BSS 


X2C 


DKH2 


FULL 


35.8447 


35.8448 


35.8439 


35.9064 


DLU(ALL)-FULL 


-0.0006 


-0.0007 


0.0040 


0.0101 


DLU(WP,NR)~FULL 


-0.1579 


-0.1579 


-0.1533 


-0.1476 


DLU(W,NR)-FULL 


0.3648 


0.3648 


0.3684 


0.3737 


DLU(WP,BP)-FULL 


1.2501 


1.2500 


1.2547 


1.2608 


DLU(W,BP)-FULL 


254.4914 


254.4914 


254.4951 


254.5018 


DLH(ALL)-FULL 


19.4202 


19.4024 


21.0255 


21.1635 


DLH(WP,NR)-FULL 


19.4229 


19.4051 


21.0284 


21.1663 


DLH(W,NR)-FULL 


19.6722 


19.6544 


21.2789 


21.4146 


DLH(WP,BP)-FULL 


20.3668 


20.3490 


21.9743 


22.1116 


DLH(W,BP)-FULL 


273.0474 


273.0296 


274.6570 


274.7922 



However, the NR approximation is much better. The errors on the reaction energy within 
the DLU(NR) scheme, which are in the range of 0.1 to 0.3 niH and thus less than 1 kJ/mol, are 
perfectly acceptable. This success is due to the correct nonrelativistic limit of the NR approxima- 
tion within the DLU scheme. Not unexpectedly, the DLH(NR) results have quite larger errors, 
because DLH(ALL) already does not provide good results for the phosphine dissociation energy 
from the WH6(PH3)3 complex. If we study the difference between DLH(NR) and DLH(ALL), we 
find that the NR approximation to actually give small errors on top of DLH(ALL). However, the 
DLH scheme is not recommended since its accuracy is not guaranteed at all. Only in those cases 
for which the DLH scheme may work, we can further add the NR approximation on top of it for 
light atoms. As can be seen from both Table IIVI and Table El different relativistic approaches 
behave basically identical with respect to the discussion of the NB and NR approximations. 
Therefore, the approximations for the reduction of computational cost discussed in this article 
can be applied to all relativistic transformation approaches. 
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VII. CONCLUSIONS 



In this work, we aimed at a rigorous local approximation to relativistic transformation 
schemes, which makes them applicable to and efficient for calculations on large molecules. We 
developed a systematic hierarchy of approximations that is based on the assembly of the unitary 
transformation from 'atomic' contributions (DLU) rather than on a local approximation directly 
applied to the matrix representation of the Hamiltonian (DLH). 

The straightforward DLH scheme, which has evolved in the development of the DKH method 
(see references given above), turns out to be not a very accurate local approximation since it only 
covers the relativistic treatment of the atom-same-atom diagonal blocks in the Hamiltonian. 
It may fail if inter-atomic distances become short so that a relativistic treatment of the off- 
diagonal 'interatomic' blocks becomes important. As a consequence, the DLH approximation 
is not suitable, for example, in studies that crucially depend on accurate electronic energy 
surfaces. Even if the molecular geometry is fixed, the DLH approach does not provide a balanced 
description for different molecular properties, since the contributions of the off-diagonal blocks 
to total observables are quite different for different properties. 

By contrast, the DLU scheme proposed by us in this paper overcomes the drawbacks of 
the DLH scheme. It is an excellent approach to take into account the atom-other-atom (i.e., 
'interatomic') off-diagonal relativistic transformations. It does not show an instability for varying 
molecular structures (i.e, for electronic energy surfaces) and properties. The size of the error 
introduced by DLU is much smaller than that of DLH; it is typically only 1% of the latter. The 
errors of the DLU approach when compared to the full relativistic transformation without any 
local approximation turn out to be tiny. They are even smaller than the difference among the 
different exact-decoupling approaches, which are already very small. 

The DLU scheme is valid for all exact-decoupling approaches, while the local X-matrix scheme 
only works for the X2C approach. Furthermore, the DLU scheme can also be applied to finite- 
order DKH approaches such as DKH2, which has been very successful in computational chemistry 
and whose computational costs are less than those of any exact-decoupling approach. 

If one has to use the DLH scheme for some reason, such as the lack of a DLU implementation 
or because of its linear-scaling behavior, the selection of a relativistic approaches is not decisive 
as the errors introduced by the DLH scheme are already higher than the difference of different 
relativistic approaches. Hence, the approximate, but fast DKH2 method may be safely combined 
with the DLH approximation, while this cannot be recommended for any exact-decoupling ap- 

26 



proach. 

However, the linear-scaling behavior of DLH is no good reason since the DLU scheme can 
also be made linear-scaling within the neighboring-atomic-blocks approximation, which produces 
negligible errors. One can further reduce the computational costs by employing the BP or NR 
approximation to all light atoms. However, the BP approximation turns out to be not suitable 
for variational calculations since it introduces large errors. Only the NR approximation yields 
acceptable results. The errors of the NR approximation will be larger if the nuclear charge of 
the atom for which the NR approximation is invoked increases. Therefore, the application of the 
NR approximation will depend on the balance of accuracy and computational cost in an actual 
calculation. 

While this work focused on the accurate approximation of exact-decoupling methods at the 
self-consistent-field level with respect to total electronic energies and first-order properties of a 
small, but well-selected test molecule set, further investigations into the DLU scheme concerning, 
for example, higher-order properties, spin-orbit and electron-correlation effects are under way 
in our laboratory. 
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